home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / tutorials / geometer / conic.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  8.6 KB  |  349 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17.  
  18. #include "parse.h"
  19. #include <math.h>
  20. #include "generic.h"
  21.  
  22. #define DIM 5
  23. #define TINY 1e-20
  24.  
  25. extern float fxmin, fxmax, fymin, fymax;
  26.  
  27. long ludecompose(double a[DIM][DIM], long indx[DIM])
  28. {
  29.     register i, j, k;
  30.     long imax;
  31.     register float sum;
  32.     float dum, big;
  33.     register double vv[DIM];
  34.  
  35.     for (i = 0; i < DIM; i++) {
  36.     big = 0.0;
  37.     for (j = 0; j < DIM; j++)
  38.         if (fabs(a[i][j]) > big)
  39.         big = fabs(a[i][j]);
  40.     if (big == 0.0) {
  41.         return 0;
  42.     }
  43.     vv[i] = 1.0/big;
  44.     }
  45.     for (j = 0; j < DIM; j++) {
  46.     if (j > 0) {
  47.         for (i = 0; i < j; i++) {
  48.         sum = a[i][j];
  49.         if (i > 0) {
  50.             for (k = 0; k < i; k++)
  51.             sum -= a[i][k]*a[k][j];
  52.             a[i][j] = sum;
  53.         }
  54.         }
  55.     }
  56.     big = 0.0;
  57.     for (i = j; i < DIM; i++) {
  58.         sum = a[i][j];
  59.         if (j > 0) {
  60.         for (k = 0; k < j; k++)
  61.             sum -= a[i][k]*a[k][j];
  62.         a[i][j] = sum;
  63.         }
  64.         dum = vv[i]*fabs(sum);
  65.         if (dum > big) {
  66.         big = dum;
  67.         imax = i;
  68.         }
  69.     }
  70.     if (j != imax) {
  71.         for (k = 0; k < DIM; k++) {
  72.         dum = a[imax][k];
  73.         a[imax][k] = a[j][k];
  74.         a[j][k] = dum;
  75.         }
  76.         vv[imax] = vv[j];
  77.     }
  78.     indx[j] = imax;
  79.     if (j != DIM - 1) {
  80.         if (a[j][j] == 0.0)
  81.         a[j][j] = TINY;
  82.         dum = 1.0/a[j][j];
  83.         for (i = j+1; i < DIM; i++)
  84.         a[i][j] = a[i][j]*dum;
  85.     }
  86.     }
  87.     if (a[DIM-1][DIM-1] == 0.0)
  88.     a[DIM-1][DIM-1] = TINY;
  89.     return 1;
  90. }
  91.  
  92. void lubksb(double a[DIM][DIM], long indx[DIM], double b[DIM])
  93. {
  94.     long i, j;
  95.     long ip, ii;
  96.     register double sum;
  97.  
  98.     ii = -1;
  99.     for (i = 0; i < DIM; i++) {
  100.     ip = indx[i];
  101.     sum = b[ip];
  102.     b[ip] = b[i];
  103.     if (ii != -1) {
  104.         for (j = ii; j < i; j++)
  105.         sum -= a[i][j]*b[j];
  106.     } else if (sum != 0.0)
  107.         ii = i;
  108.     b[i] = sum;
  109.     }
  110.     for (i = DIM-1; i >= 0; i--) {
  111.     sum = b[i];
  112.     if (i < DIM-1) {
  113.         for (j = i+1; j < DIM; j++)
  114.         sum = sum - a[i][j]*b[j];
  115.     }
  116.     b[i] = sum/a[i][i];
  117.     }
  118. }
  119.  
  120. void linetangenttoconic(line *l)
  121. {
  122.     conic *c = (conic *)l->c.p1;
  123.     vertex *v = (vertex *)l->c.p2;
  124.     if (v->w == 0.0) {
  125.         l->v1.xw = l->v2.xw = 1000.0;
  126.     l->v1.yw = -1.0; l->v2.yw = 1.0;
  127.     l->v1.w = l->v2.w = 0.0;
  128.     homogenizeline(l);
  129.     return;
  130.     }
  131.     if ((v->c.type == VERTONCONIC && (v->c.p1 == c)) ||
  132.         (v->c.type == VERTLINECONIC1 && (v->c.p2 == c)) ||
  133.     (v->c.type == VERTLINECONIC2 && (v->c.p2 == c))) {
  134.     double x = v->xw, y = v->yw;
  135.     // v had better be on the conic
  136.     l->v1.xw = x; l->v1.yw = y; l->v1.w = 1.0; l->v2.w = 1.0;
  137.     float dx = (c->bb*x+2.0*c->cc*y+c->ee);
  138.     float dy = - (c->bb*y+2.0*c->aa*x+c->dd);
  139.     float den = fabs(dx) + fabs(dy);
  140.     l->v2.xw = x + dx/den;
  141.     l->v2.yw = y + dy/den;
  142.     homogenizeline(l);
  143.     } else {
  144.     l->v1.w = l->v2.w = 1.0;
  145.     float alpha = c->aa*v->xw + c->bb*v->yw/2.0 + c->dd/2.0;
  146.     float beta = c->bb*v->xw/2.0 + c->cc*v->yw + c->ee/2.0;
  147.     float gamma = c->dd*v->xw/2.0 + c->ee*v->yw/2.0 + c->ff;
  148.     float A = c->aa*beta*beta - c->bb*alpha*beta + c->cc*alpha*alpha;
  149.     float B = 2*c->cc*alpha*gamma - c->bb*beta*gamma + c->dd*beta*beta - c->ee*alpha*beta;
  150.     float C = c->cc*gamma*gamma - c->ee*gamma*beta + c->ff*beta*beta;
  151.     float Disc = (B*B - 4*A*C);
  152.     if (Disc < 0) {
  153.         l->v1.xw = l->v2.xw = 1000.0;
  154.         l->v1.yw = -1.0; l->v2.yw = 1.0;
  155.         l->v1.w = l->v2.w = 0.0;
  156.         homogenizeline(l);
  157.         return;
  158.     }
  159.     Disc = sqrt(Disc);
  160.     if (l->c.type == LINETANGENTCONIC1)
  161.         l->v1.xw = (-B + Disc)/(2.0*A);
  162.     else
  163.         l->v1.xw = (-B - Disc)/(2.0*A);
  164.     l->v1.yw = (-gamma - alpha*l->v1.xw)/beta;
  165.     l->v2.xw = v->xw; l->v2.yw = v->yw;
  166.     homogenizeline(l);
  167.     }
  168. }
  169.  
  170. void conic_vvvvv(conic *c)
  171. {
  172.     double x[5], y[5], b[5];
  173.     double a[5][5];
  174.     long indx[5];
  175.     
  176.     vertex *v1 = (vertex *)c->c.p1;
  177.     vertex *v2 = (vertex *)c->c.p2;
  178.     vertex *v3 = (vertex *)c->c.p3;
  179.     vertex *v4 = (vertex *)c->c.p4;
  180.     vertex *v5 = (vertex *)c->c.p5;
  181.     x[0] = v1->xw; y[0] = v1->yw;
  182.     x[1] = v2->xw; y[1] = v2->yw;
  183.     x[2] = v3->xw; y[2] = v3->yw;
  184.     x[3] = v4->xw; y[3] = v4->yw;
  185.     x[4] = v5->xw; y[4] = v5->yw;
  186.  
  187.     for (long i = 0; i < 5; i++) {
  188.     a[i][0] = x[i]*x[i];
  189.     a[i][1] = x[i]*y[i];
  190.     a[i][2] = y[i]*y[i];
  191.     a[i][3] = x[i];
  192.     a[i][4] = y[i];
  193.     }
  194.     b[0] = b[1] = b[2] = b[3] = b[4] = -1.0;
  195.         
  196.     if (0 == ludecompose(a, indx)) {
  197.     c->aa = c->bb = c->cc = c->dd = c->ee = c->ff = 0.0;
  198.     return;
  199.     }
  200.     lubksb(a, indx, b);
  201.  
  202.     c->aa = b[0];
  203.     c->bb = b[1];
  204.     c->cc = b[2];
  205.     c->dd = b[3];
  206.     c->ee = b[4];
  207.     c->ff = 1.0;
  208. }
  209.  
  210. void conic_lllll(conic *c)
  211. {
  212.     double x[5], y[5], w[5], b[5];
  213.     double a[5][5];
  214.     long indx[5];
  215.     
  216.     line *l1 = (line *)c->c.p1;
  217.     line *l2 = (line *)c->c.p2;
  218.     line *l3 = (line *)c->c.p3;
  219.     line *l4 = (line *)c->c.p4;
  220.     line *l5 = (line *)c->c.p5;
  221.     x[0] = l1->XW; y[0] = l1->YW; w[0] = l1->W;
  222.     x[1] = l2->XW; y[1] = l2->YW; w[1] = l2->W;
  223.     x[2] = l3->XW; y[2] = l3->YW; w[2] = l3->W;
  224.     x[3] = l4->XW; y[3] = l4->YW; w[3] = l4->W;
  225.     x[4] = l5->XW; y[4] = l5->YW; w[4] = l5->W;
  226.  
  227.     for (long i = 0; i < 5; i++) {
  228.     a[i][0] = x[i]*x[i];
  229.     a[i][1] = x[i]*y[i];
  230.     a[i][2] = y[i]*y[i];
  231.     a[i][3] = x[i]*w[i];
  232.     a[i][4] = y[i]*w[i];
  233.     b[i] = -w[i]*w[i];
  234.     }
  235.     //b[0] = b[1] = b[2] = b[3] = b[4] = -1.0;
  236.         
  237.     if (0 == ludecompose(a, indx)) {
  238.     c->aa = c->bb = c->cc = c->dd = c->ee = c->ff = 0.0;
  239.     return;
  240.     }
  241.     lubksb(a, indx, b);
  242.     //float det = b[0]*b[2] + b[1]*b[3]*b[4]/4 - b[2]*b[3]*b[3]/4
  243.     //            - b[0]*b[4]*b[4]/4 - b[1]*b[1];
  244.     c->aa = (b[2] - b[4]*b[4]/4);
  245.     c->bb = (b[3]*b[4]/2 - b[1]);
  246.     c->cc = (b[0] - b[3]*b[3]/4);
  247.     c->dd = (b[1]*b[4]/2 - b[2]*b[3]);
  248.     c->ee = (b[1]*b[3]/2 - b[0]*b[4]);
  249.     c->ff = (b[0]*b[2] - b[1]*b[1]/4);
  250.     //c->aa = b[0];
  251.     //c->bb = b[1];
  252.     //c->cc = b[2];
  253.     //c->dd = b[3];
  254.     //c->ee = b[4];
  255.     //c->ff = 1.0;
  256. }
  257.  
  258. void projecttoconic(vertex *v, conic *con)
  259. {
  260.     float x = v->xw, y = v->yw;
  261.     float dist = 1.0e30, dt;
  262.     float xval, yval, xbest, ybest, u;
  263.     
  264.     for (xval = fxmin; xval <= fxmax; xval += .01) {
  265.     float a = con->cc;
  266.     float b = con->bb*xval + con->ee;
  267.     float c = con->aa*xval*xval + con->dd*xval + con->ff;
  268.     if ((u = b*b - 4.0*a*c) >= 0) {
  269.         float disc = sqrt(b*b - 4.0*a*c);
  270.         float y1 = (-b + disc)/(2*a);
  271.         float y2 = (-b - disc)/(2*a);
  272.         if ((dt=(x-xval)*(x-xval)+(y1-y)*(y1-y))<dist) {
  273.             dist = dt;
  274.         xbest = xval;
  275.         ybest = y1;
  276.         }
  277.         if ((dt=(x-xval)*(x-xval)+(y2-y)*(y2-y))<dist) {
  278.             dist = dt;
  279.         xbest = xval;
  280.         ybest = y2;
  281.         }
  282.     }
  283.     }
  284.     for (yval = fymin; yval <= fymax; yval += .01) {
  285.     float a = con->aa;
  286.     float b = con->bb*yval + con->dd;
  287.     float c = con->cc*yval*yval + con->ee*yval + con->ff;
  288.     if ((u = b*b - 4.0*a*c) >= 0) {
  289.         float disc = sqrt(b*b - 4.0*a*c);
  290.         float x1 = (-b + disc)/(2*a);
  291.         float x2 = (-b - disc)/(2*a);
  292.         if ((dt=(x-x1)*(x-x1)+(y-yval)*(y-yval))<dist) {
  293.             dist = dt;
  294.         xbest = x1;
  295.         ybest = yval;
  296.         }
  297.         if ((dt=(x-x2)*(x-x2)+(y-yval)*(y-yval))<dist)  {
  298.             dist = dt;
  299.         xbest = x2;
  300.         ybest = yval;
  301.         }
  302.     }
  303.     }
  304.     for (xval = xbest - .01; xval <= xbest + .01; xval += .0003) {
  305.     float a = con->cc;
  306.     float b = con->bb*xval + con->ee;
  307.     float c = con->aa*xval*xval + con->dd*xval + con->ff;
  308.     if ((u = b*b - 4.0*a*c) >= 0) {
  309.         float disc = sqrt(b*b - 4.0*a*c);
  310.         float y1 = (-b + disc)/(2*a);
  311.         float y2 = (-b - disc)/(2*a);
  312.         if ((dt=(x-xval)*(x-xval)+(y1-y)*(y1-y))<dist) {
  313.             dist = dt;
  314.         xbest = xval;
  315.         ybest = y1;
  316.         }
  317.         if ((dt=(x-xval)*(x-xval)+(y2-y)*(y2-y))<dist) {
  318.             dist = dt;
  319.         xbest = xval;
  320.         ybest = y2;
  321.         }
  322.     }
  323.     }
  324.     for (yval = ybest - .01; yval <= ybest+.01; yval += .0003) {
  325.     float a = con->aa;
  326.     float b = con->bb*yval + con->dd;
  327.     float c = con->cc*yval*yval + con->ee*yval + con->ff;
  328.     if ((u = b*b - 4.0*a*c) >= 0) {
  329.         float disc = sqrt(b*b - 4.0*a*c);
  330.         float x1 = (-b + disc)/(2*a);
  331.         float x2 = (-b - disc)/(2*a);
  332.         if ((dt=(x-x1)*(x-x1)+(y-yval)*(y-yval))<dist) {
  333.             dist = dt;
  334.         xbest = x1;
  335.         ybest = yval;
  336.         }
  337.         if ((dt=(x-x2)*(x-x2)+(y-yval)*(y-yval))<dist)  {
  338.             dist = dt;
  339.         xbest = x2;
  340.         ybest = yval;
  341.         }
  342.     }
  343.     }
  344.     if (dist < 10000000.0) {
  345.     v->xw = xbest; v->yw = ybest; v->w = 1.0;
  346.     } else v->w = 0.0;
  347. }
  348.  
  349.